#include "../Random.hpp"
#include <armadillo>
#include <boost/math/distributions/normal.hpp>
#include "../utils.hpp"


#define _M 100



class Probit_MS : public Density::GeomBridge
{
		friend class Density::GeomBridge;
	public:
		Probit_MS(double *Y, mat X, Distribution::Distribution *S,mat s) : Density::GeomBridge(Y, X, S){ 
			_s=s;
			Set_p(X.n_cols-1);
			_is=inv(s);
			fL=new Laplace;
		}
		~Probit_MS(){
			delete fL;
		}
		inline double Likelihood(mat gamma)
		{
			int p=Get_p();
			mat X=Get_X();
			mat gammat(p+1,1);
			mat h=sum(gamma,0);
			double res=0;
			int nb=h(0,0)+1;
			int n=X.n_rows;	
			mat X2(n,nb);
			cout << "test"<< p+1;
			int kk=1;
			for(int i=0;i<(p+1);i++)
			{
					
				if(i!=0){
					gammat(i,0)=gamma(i-1,0);
					if(gamma(i-1,0)!=0)
					{
						X2(span::all,kk)=X(span::all,i);
						kk++;
					}

				}else{
					gammat(i,0)=1;
					X2(span::all,0)=X(span::all,0);
				}
			}
		
			mat Xt=X*diagmat(gammat);
			mat m(nb,1);
			m.zeros();
			double *Y=Get_Y();
			mat foo(nb,nb);
			foo.fill(10);
			mat s=diagmat(foo);
		//	cout << X2;	
			(*fL)(X2,Y,m,_s);			
			mat Sig=(*fL).Get_Sig();	
			mat Mu=(*fL).Get_Mu();
			int k=0;
			int g=0;
			mat m1(1,p+1);
			m1.fill(0);
			
			mat gammaS=gammat*gammat.t();
			/*mat S=randu<mat>(p+1,p+1);
				
			for(int i=0;i<p+1;i++)
			{
				if(gammat(i,0)!=0)
				{
					m1(0,i)=Mu(k,0);
					g=0;
					for(int j=0;j<p+1;j++)
					{
						if(gammat(j,0)!=0)
						{
							S(i,j)=Sig(k,g);
							g++;
						}			
					}
					k++;				
				}	
			}
			cout << S;
		*/
			//Distribution::Gaussian G(Xt.n_cols,m,s);
			Distribution::Gaussian G(nb,Mu,Sig);
			
			mat prop=G.r();
			//cout << "test" << nb-prop.n_rows;
			double sum=0;
			for( int i=0;i<_M;i++)
			{
				prop=G.r();
				mat beta=add_mat(prop,gammat,p+1).t();
				sum+=exp(Lik_Prob(beta,Xt)+Priorb(beta)-G.d(prop));
			}
			res=log(sum/_M);
			//cout << res << "\n";
			

			return res;

		}
		double Lik_Prob(mat theta1,mat Xt)
		{
			mat theta=theta1.t();
			double L=0;
			double sum=0;
			int m=theta.n_cols;
			//Xb
			int n=Get_n();
			for(int i=0;i<n;i++)
			{	
			//	cout << theta.col(0) << "\\";
				mat foo= Xt.row(i);
				sum=0;
				for(int j=0;j<m;j++)
				{
					sum+=theta(0,j)*foo(0,j);
				}
				if(Get_Y()[i]==1){
					L+=log(Phi(sum));
				}else if(Get_Y()[i]==0){
					L+=log(1-Phi(sum));
				}else{
					cout << "Warning values of Y!= {0,1}";
				}
			}
			//cout << exp(L);
			return L;

		}
		double Priorb(mat theta)
		{
			double res=0;
			mat foo=theta.t();
			mat bar=foo*_is*foo.t();
			//cout << "//" << bar;
			res=-0.5*log(2*PI*det(_s))-0.5*bar(0,0);
			return res;

		}

		double Prior(mat theta){
			return -0.693;
		}

	private:
		Laplace *fL;
		mat _s;
		mat _is;
	
};
